Multiple-Point and Multiple-Time Correlations Functions in a Hard-Sphere Fluid 
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A recent mode coupling theory of higher-order correlation functions is tested on a simple hard- 
sphere fluid system at intermediate densities. Multi-point and multi-time correlation functions of 
the densities of conserved variables are calculated in the hydrodynamic limit and compared to 
results obtained from event-based molecular dynamics simulations. It is demonstrated that the 
mode coupling theory results are in excellent agreement with the simulation results provided that 
dissipative couplings are included in the vertices appearing in the theory. In contrast, simplified mode 
coupling theories in which the densities obey Gaussian statistics neglect important contributions to 
both the multi-point and multi-time correlation functions on all time scales. 
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I. INTRODUCTION 



Over the last few years, the emergence of multi- 
dimensional NMR and non-resonant non-Hnear Ra- 
man [^|-^ techniques has generated renewed interest in 
the information content of higher-order correlation func- 
tions involving time-correlations of dynamical quantities 
at multiple points and time separations. These exper- 
imental developments hold great promise for the eluci- 
dation of the nature of the underlying dynamics giving 
rise to complex relaxation behavior in super-cooled liq- 
uids, polymeric systems, and proteins |^-^. Concur- 
rently, simulation studies probing the microscopic ori- 
gin of dynamical heterogeneity in dense systems have 
made use of the increased information content available 
in multiple-point [ p^ and multiple-time JTTI] correlation 
functions. 

Although there has been some recent work attempt- 
ing to reproduce simulation results for the off-resonant 
fifth-order Raman response function JT2|-p^, there has 
been little theoretical work to establish a microscopic 
theory for general higher-order correlation functions. In 
a previous article |15| , a general mode coupling theory 
was developed in which the long time behavior for multi- 
point and multi-time correlation functions was expressed 
in terms of ordinary two-time, two-point correlation func- 
tions of a set of slow variables which are coupled by ver- 
tices containing both static (called Euler) and dynamic 
(called dissipative) correlations. The theory is based 
upon the assumption that the long-time dynamics of ar- 
bitrary variables is a functional of a set of slow modes 
of the system. The long-time dynamics of higher-order 
correlation functions is then described by isolating the 
component of the relevant variables along multi-linear 
products of the slow variables, resulting in expressions for 
the higher-order correlation functions in terms of the sum 
of an infinite number of multi-point correlation functions 
of slow modes. The formulation is made tractable by a 
cumulant expansion method (called iV-ordering [|l6|,0) 
in which multi-point correlation functions are factored 



into convolutions over the familiar two-point, two-time 
correlation functions of the slow modes. In this way, 
the need to simplify the mode-coupling expressions for 
higher-order correlation functions based on an assump- 
tion of Gaussian statistical behavior of the slow modes is 
avoided. It was suggested that simple mode coupling the- 
ories [|l8|,^ based upon this Gaussian assumption lead 
to a relatively poor description of the long time behavior 
of higher-order correlation functions. 

The purpose of this article is to validate the mode 
coupling theory expressions for multi-point and multi- 
time correlation functions by examining the simplest non- 
trivial system, the hard-sphere liquid. The hard-sphere 
liquid is a very useful system to examine theoretically 
since the simple form of the interaction potential allows 
static correlation functions to be related to the radial dis- 
tribution function at contact. In turn, the radial distri- 
bution function can be approximated using an accurate 
equation of state, such as that of Carnahan and Starling 
POI , which relates the pressure to the density and the 
temperature. In addition, excellent predictions exist for 
dynamical properties of hard sphere systems based on de- 
tailed kinetic theory . Another advantage of looking 
at hard sphere systems is that the dynamics of the sys- 
tem can be simulated very efficiently using event-based 
molecular dynamics methods since particles evolve 
freely between collisions, thereby allowing good statistics 
to be obtained from simulations. 

We shall focus on systems of moderate reduced densi- 
ties (p* = 0.25) in which "mode coupling" effects leading 
to non-exponential relaxation of correlation functions of 
linear densities, such as the dynamical structure factor, 
can be neglected. In particular, we target correlation 
functions of long-wavelength fluctuations which decay on 
long time scales and exhibit complicated higher-order 
correlation functions. 

This paper is organized as follows: In Section II, the 
mode coupling formalism developed in Ref. jist is re- 
viewed and adapted to the hard sphere system. Explicit 
expressions are presented for three-point and three-time 
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correlation functions of involving linear densities of num- 
ber (or mass), transverse, and longitudinal velocities. In 
Section III, simulation methods particularly suited for 
calculating higher-order correlation functions in a hard- 
sphere system are discussed. In Section IV, the predic- 
tions of the mode coupling theory are compared to the 
simulation results for relatively simple three-point and 
three-time correlation functions, and it is demonstrated 
that dissipative parts of vertices provide additional im- 
portant couplings to those at Euler order. The results 
are contrasted with those obtained within the framework 
of Gaussian mode coupling theory [|l8|,^ . Finally, con- 
clusions of the study are given in Section V. 



II. THEORETICAL FORMULATION 



The system under consideration is composed of TV 
particles of mass m and diameter a in a volume V = 
Lz- The particles interact through the two- 
body hard-sphere potential, 

y(r) = |o (1) 

^ ' I oo II r < a. ^ ' 

Given the form of the potential, the dynamics generated 
by the Hamiltonian conserves the total number of parti- 
cles iV, the total angular momentum, the linear momenta 
P, and the energy E of the system. In Ref. |l^], expres- 
sions for the long time behavior of correlation functions 
were obtained under the assumption that the slowly vary- 
ing part of an arbitrary dynamical variable is an analytic 
function of a set of slow variables A of the system. An 
essential part of successfully applying the formalism to a 
particular system is the identification of a complete set 
of slow variables. To identify the slow modes of the sys- 
tem, it is helpful to consider the local densities of the 
conserved variables iV, P and E, 

N 

7V(r) = ^J(r-r,), 

i=i 

N 

P(r) = ^p,5(r-r,), 

i=l 

W / 2 -1 \ 

where and Pi are the spatial position and momentum 
of particle i. Noting that the Fourier transform of these 
densities. 



N 

i=l 

N 

Pk-^p.e*-', (2) 

i=l 

are slowly varying quantities for small fc — |k| since their 
time derivatives are proportional to fc, the minimal set of 
slow variables must include all the "hydrodynamic" 
variables {A'k, Pk, £'k} with fc smaller than some cut-off 
wave vector fc^. For our purposes, it is convenient to 
work with a slightly different basis set Ak, composed of 
the variables A^k, = P£, Ty^ = P^, Tsk = Pk, and 
ifk = (SA^k - 2/3£:k)/\/6, where /? = l/(fcBr) is the in- 
verse temperature of the system, P^ is the x-component 
of the vector Pk, and k is taken along the x-axis. Note 
that the Tik and T2V. are the transverse modes of the mo- 
mentum density, while Lk is the longitudinal momentum 
density. With this definition of basis set, the matrix 

is diagonal in the hydrodynamic labels a and 6, where 
(• • • ) denotes the grand-canonical ensemble average. The 
non-linear dependence of the dynamical variables is ex- 
pressed in terms of a "multi-linear" basis set 

Qo = 1 

Qi = Ak - (Ak> = ik (3) 

Q2 = Qk-qQq — (Qk-qQq) — (Qk-qQqQ*) ' ' Ql, 



where the "•" notation denotes a sum over components of 
the column vector ^k (the indices of the hydrodynamic 
variables iVk, Lk, Tik, T2k, and The subtractions 

in the basis set defined in Eq. (||) are included to ensure 
that the multi-linear matrix, 

Ki^ = (QiQ:„) = {QiQ*J Stm, (4) 

is diagonal in mode-order I. The slow part of any dynam- 
ical variable C is removed by the projection operator 

00 

vc ^J2(^Qi)^a'Qi^ (5) 

and the complementary projection operator = 1 — V 
projects onto the orthogonal sub-space. 

Writing the three-point correlation function 
(^k-q(i)^q(^)^-k) in terms of the basis set, we obtain 

(ik-q(<)iq(<)i-k) = (ik(t)i-k) • K^l ■ (i-kik-qiq) 
+ G^-q,q;k(i) ' (6) 
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where G"""(t) 



^k-q,q;k(^) 



(Qm(t)Q 

(Q2(k 



* -ft'n,! , and, in particular, 



q, q;i)^-k> 



(7) 



Note that Eq. (g) is exact in the Umit t — > by con- 
struction of the basis set. Utihzing projection operator 
techniques 
the mu: 



4yi7| ] and cumulant expansion methods 
ti-point correlation function G^_^ q k(^) '^^^^ 
be expressed in terms of two-point, two-time correlation 
hmctions as lIlSl, 



^k-q,q;k(0 



f Gi\{t-T)G]^{t 
Jo 



r) 



^'^k-q.q;k 



Gii(r)dr, 



(8) 



where G^{t) = (ylk(i)^-k)/(^k^-k) are the normal- 
ized, two-point and two-time correlation functions of the 
linear densities, and the "vertices" are given by 







K- 



with the fluctuating force 0; (t) defined by 



(9) 



(10) 



where C is the Liouville operator. 

Similarly, it can be shown that the three-time correla- 
tion function. 



G"i(ti,t2) = {Qiih + t2)Qiih)Ql) ■ 
can be approximately written as p^ ], 

G"i(t2,ti) =G"(i2)*Mi"*G"(ti) 
+Gi2(t2)*M2ii*G"(ti) 



+G''it2) * M''' * G'^ih) + 0{N-'), 
where M'™" is given by 



K- 



(11) 



(12) 



(13) 



Furthermore, it was shown in Ref. |15| that G {12) can 



be written in terms of the two-point, two-time functions 
and the AP^ vertices in a manner analogous to Eq. (^). 

The symmetry properties of the Hamiltonian can be 
used to greatly simplify the analysis of higher-order cor- 
relation functions. For example, since the Hamiltonian 
Ti. is invariant under the transformation TTi. = H, where 
the self- adjoint time-reversal operator T acts on an arbi- 
trary phase point (r^,p^) by T(r^,p^) = (r^,— p^), 
all time correlation functions considered here have well- 
defined symmetry properties under T: TA^ = 7a^ki 
where 7a = 1 for a ~ N, H and 7a = — 1 for a = 
Ti,T2,L. Furthermore, since the Liouville operator C 
transforms as TC = —CT, it is easy to show [ p"5|j25| 

that (^£(i)^^_k> = 7a7fc(^£(-*)^-k>- I* is straight- 
forward to extend these arguments to multi-time corre- 
lation functions for which {Al_^{ti + t2)^q(ti)^lk) = 

lalblc{AU-tl-t2)A'^{-tM-k)- 



A. Three-point correlations 



We now turn our attention to evaluating the expres- 
sions for three-point correlation functions of three lin- 
ear densities of the form in Eq. (^) for several differ- 
ent combinations of wave-vector and hydrodynamic la- 
bels in terms of the linear-linear correlation functions 
G^^{t). For simplicity, we consider correlation functions 
involving the transverse momentum mode r2k henceforth 
abbreviated as just Tk- From symmetry considerations 
it is easy to establish that the linear-linear correlation 
function G^''(t) = unless a = T, which simplifies the 
subsequent analysis. 

Looking first at the correlation function, 
(rk-q(t)Tq(i)iV_k), using Eq. (I) we have. 



(rk-q(t)rq(t)7V-k) 



{mt)A^u 



-(i^k^k-qTq) 



+ Gj^q~q;k(0(^kiV_k) 



(14) 



where the repeated index a is summed over the labels for 
iV, T, L and H, and 



^k-q,q; 



u^- (Qr(k-q,q,t)iV,k) 



The replacement of the "21" super-script in G^^^ ^.^(t) 
by "TT; A^" above is meant to denote the specific hy- 
drodynamic labels under consideration. The semi-colon 
separating the labels indicates that the labels "TT" cor- 
respond to the bi-linear density, whereas the "A'^" labels 
the linear density in Eq. (0). Noting that (i^klk-qTq) 
vanishes unless a — N,H, the first part of Eq. ( p^ can 
be written as. 



S{k) 



-(7V_kTk-qTq/ 



(iVk(t)g_k) 

(N) 



(-H^-k?k-qTq) 



= -(7Vk(t)iV_k)-^(iVk(0^-k), 



where S{k) = (iVkiV_k) is the static structure fac- 
tor. The normalized multi-point correlation function 
Gk-q^ k(^) basis set variable Q'^'^ can be eval- 

uated using Eq. (||), 

Grq^q;kW= fGr_^{t-r)Gl^{t~r) 

JO 

Mn^\MGt''{r)dr, (15) 

where a is summed over the labels N and H only since 
the M'^'^''°- vertex vanishes when a = LorT. The explicit 
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form of the vertix is given by Eq. (H), which involves an 
"Euler part" , 



ATT Aa 
•^k-q.q^-k 



{AlA-_^) 

and a dissipative part 

) 

These contributions to the vertex can be evaluated as de- 
tailed in the appendix, and one finds that to leading order 
in the wave-vector, only the a = L term contributes at 
order k (with corrections of order k^), whereas the other 
vertex with a = H gives a contribution proportional to 
k^. Putting all this together, we obtain the expression, 

(rk-q(t)rq(i)A^-k> = ^(A^kWiv_o - -^{N^m-^) 

(16) 



where 



'^k-q,q;ki*-' 



and the functions Gj'^(r), G^^(t), and G^^(r) are 
given explicitly by 

Crir) - {T^{r)T^^)/m{N)kBT 
G^M - (Lk(r)iV_k>/^(fc) 
(r) - (ffk(r)iV_k>/^(fc)- 



The vertices M^^^^.^^ and M^^'^ 



''k-q^q;k given in the 
appendix. Note that if the dissipative parts of the ver- 
tices are ne glec ted, only the first time convolution inte- 
gral in Eq. (0) contributes to (Tk_q(t)Tq(t)iV_k)- 

From similar considerations, it is not difficult to obtain 
expressions for other correlation functions. For example, 
we find the multi-point functions (Tk-q(t)Lq(t)r_k) and 
(rk-q(t)iVq(t)T_k) are given by 



(rk-q(t)iq(t)T_k) = G^-q',q;k(0 



(18) 



(rk-q(t)iVq(t)T_k) = ^(rk(t)r_k) +G^i^'q,kW (19) 



(N) 



with 



" / ^k-q 

Jo 



\{t-r) 



(20) 



E G-(t-r)M-:^^^,Gr(r)dr, 



and 



GL^^;kW= / G-q(t-r) 



(21) 



G-%t-r)Ml-^^^^Gr{r)dr 



a=L,N,H 



In Eqs. (|2^) and ( ^l|) , to leading order in the wave- 
vectors, the vertex M^^'^ for a — L contributes at Eu- 
ler order and is imaginary, while the vertices M'^^''^ and 
]\^TN;T contribute at dissipative order and are real. How- 
ever, since G^^{t - r), G^^{t - r), G^^{t - r) and 
Gq"{t - t) are real and G^^(i - r), G^^{t ~ r) and 
Gq^ (t—T) are purely imaginary by time-reversal symme- 



try, the correlation functions G 



are real whereas G^f'^q.jj.(t) is purely imaginary. Note 
that at t = Q, the expressions for the functions G^^ in 
Eqs. (p7|), (|2l| ) and ( po|) vanish and the multi-point cor- 
relation functions are given exactly. 



TT-N 
k— q,q;k 



[t) and G 



TN;T 
k— q,q;k 



it) 



B. Three-time correlations 

The three-time correlation functions, 

(Tk-q(il +t2)iq(tl)r_k) 



G'k-q,q,k(*l' ^2) — 



GZ-q.q.kitl^h) 



(N) rnksT 

(rk_q(tl +t2)iVq(<l)r_k) 



(N) rnksT 



can be evaluated in a straightforward fashion using the 
results of the previous section. In Ref . ||l^ , it was shown 
that the multi-time vertices m'™" reduce to very sim- 
ple forms to leading TV-order, with corrections of order 
M/N ^ kcU « 10~^ for systems of moderate density 
p^t . Using the reduced forms of M^-^^ and M^^^ and 
Eq. (jlj) , the leading A^-order expressions for these multi- 
time functions are. 



Gk^q,q,k(*li*2) - G^.'^^^ _^{t2) (LqL^q) G(^^ {tl 



Gi-q{t2) Gkf'q^q:k(*l)> 



and 



Gk-q,q,kitl^t2) 



G'k^q(i2 



{N) 



(22) 



(23) 



+ G^5q^k,-q(^2) S{q) Grit,) 
+ G'k-q(^2) Gk^qlq.k(il). 

Using the symmetry properties of G''™'{t) ||l^, one 
write G'^;^'^(t)(»*) = {G^^'-^i-t))*, and Eqs. d 
and (p3) can be expressed in terms of G^^(t) alone as. 



G 



TLT 



k-q,q,k(*l:*2) — -G!^kiq;q-k(^2) Gk (^l) 

TL;T 



+Gk-q{t2) Gk_q_q.k 



ih), (24) 
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and 



Gk^<^^q,k(*l: ^2) 



or 



Siq) 



{N) 



(25) 



+ Gk-q(i2) Gk-q,q;k(^l)> 

where the time-reversal symmetry properties and the be- 
havior under complex conjugation of the G^^{t) correla- 
tion functions has been used. 

In Section IV, these expressions will be compared to 
results from simulations of a hard sphere system at mod- 
erate density. 



III. SIMULATION METHOD 

The dynamics of hard spheres consist of free, rectilin- 
ear motion until the distance between two spheres (i and 
j) becomes equal to their diameter a, at which point an 
instantaneous collision takes place, leading to the mo- 
mentum changes 

p, Pj - (T[(pi - Pj) • a] 
Vj Pj + o-[(Pj: - Pj) • <t], 

where the collision normal a equals (r j — ) /a at contact. 

Due to the simplicity of the equations of motion, the 
dynamical evolution of the hard sphere system can be 
computed exactly using an event-driven procedure in 
which one calculates the first possible collision of all 
spheres under the assumption that no other particles col- 
lide. The phase point of the system is then evolved up to 
the time of the earliest of these collisions, and the process 
is repeated until the total desired run time is completed. 

Without additional bookkeeping, the number of 
spheres with which a specific particle can collide is A^ — 1, 
and hence 0{N) calculations of collision times are re- 
quired for each particle after it collides. As the number of 
collisions per unit time is extensive, the simulation time 
scale increases quadratically with the number of parti- 
cles. Considerable improvements in simulation efficiency 
can be gained using a division of the system into regions 
(called cells) and data structures to optimize the search 
for the next collision time 

To use the cell structure in a simulation, the system of 
dimension x Ly x Lz is divided into an integer num- 
ber of cells of dimension Ix x ly x 1^, where each of the 
lengths Ix, ly and 1^ is no smaller than the diameter of 
the hard spheres. Now, in addition to the collision events 
between spheres, the cell in which each sphere is located 
and the time at which the particle will leave its cell is 
recorded. This is advantageous because the number of 
spheres that can collide with a given sphere before a par- 
ticle moves out of its cell is proportional to the number 
of spheres in its vicinity, i.e., the spheres in the same 
cell or in one of the 26 neighboring cells. Using the cell 



structures, the number of spheres within the vicinity of 
a given particle is of order 0{lxlylzN / LxLyL^) = 0(1), 
provided the lengths of each cell are on the order of the 
diameter of the particles, and hence far fewer collision 
times of pairs of particles must be computed after each 
collision event. However, the use of cells comes at the 
cost of increasing the complexity of the event driven sim- 
ulation since after a crossing event for a given particle, 
the collision times of the given sphere with spheres that 
previously were not in its vicinity must be considered, 
and, if necessary, the first stored collision event adjusted. 
In addition, the next crossing time in the same direction 
is re-calculated. Similarly, after a collision event between 
two particles, new collisions within the same cell as well 
as the new cell-crossing times must be calculated for the 
particles involved in the event. 

Even though the calculations after a crossing or a col- 
lision event are of 0(1) when many cells are used, it is 
still necessary to search the event list of each sphere to 
find the earliest event in the simulation. If the spheres 
are simply stored in a linear array, this implies a look-up 
time that scales linearly with N , and the algorithm scales 
as N"^ as before, though with a considerably lower pref ac- 
tor than without the cell structure. If, on the other hand, 
the spheres are stored in a binary tree, ordered according 
to their first event, the search for the first event scales 
as the logarithm of the number of elements, which in our 
implementation is N . Deleting an element from the tree 
is an 0(1) operation, while the insertion of new elements 
into the tree requires a tree search, which scales as InA^. 
Since the number of crossing and collision events is ex- 
tensive, the algorithm scales as A^lniV, and the overall 
speed up of the algorithm over a simple event-based sim- 
ulation behaves as N/ In TV. It should be noted, however, 
that the cell structure reduces the number of collisions to 
be considered to a large extent, so the prefactor is also 
quite reduced. 

There is some flexibility in selecting the size of the 
cells to be used in the simulation. Larger cells require 
fewer crossing times to be calculated at the expense of 
increasing the number of collisions which must be com- 
puted within each cell. As Rapaport has noted js^, the 
optimal choice of the dimensions of the cell for systems 
of low density is intermediate between the size of the full 
system and the diameter of the hard spheres, whereas 
the smallest possible cells make for the fastest simula- 
tion for higher densities. In the simulations reported in 
the next section, the optimal length of cells was found to 
correspond roughly with the diameter of the hard sphere 
particles. 



IV. RESULTS AND DISCUSSION 

In this section, the mode coupling expressions for the 
higher-order correlations functions given in section II are 
compared to those obtained from event-based molecu- 
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lar dynamics simulations in the micro-canonical ensem- 
ble at an inverse temperature (3 — 3. The size of 
the periodic system in the simulation was chosen to be 
= Ly = = 15.7526, such that for N = 1382 
hard sphere particles of diameter a = 1, the reduced 
density p* = 0.25 {p/pc where pc is the density at close- 
packing) and the magnitude of the smallest wave-vector 
koa = 2Tia/Lx — 0.398867 coincide with one of the cases 
in Ref. |31|. A total number of 15"^ = 3375 cells were 
used, leading to a collision rate (including data collec- 
tion) of roughly 3.2 x lOVhour of CPU time on a 600 
MHz digital 21164 processor. The event-dynamics sim- 
ulations were run on a 9 nodes of a 30-node "Beowulf" 
cluster for a total of 4, 100 CPU hours, where each node 
carried out 3, 750 short molecular dynamics trajectories 
of approximately 403, 500 collisions. The initial configu- 
ration of the system for each of the individual runs was 
randomly chosen using a simple rejection method. In all 
results reported below, time is expressed in dimensionlcss 
units t/tc, where tc is the mean collision time calculated 
from the simulation. At the density and temperature of 
the simulation, the mean collision time is roughly half 
the time tm it takes a particle to move over a distance 
equal to its diameter a {tc/tm ~ 0.412). 

To evaluate time correlation functions in the simula- 
tion, the values of the linear densities Ak(t) were calcu- 
lated for a set of wave- vectors at M+1 fixed time intervals 
t — 0, At, 2 At, . . . and stored in an array where 
the index k runs over the wave vector indices and i runs 
from to M. In all molecular dynamics runs, the time in- 
terval At/tra = 0.15, and M = 400. Two-point, two-time 
correlation functions for a given time interval were accu- 
mulated on the fly by storing the product of accumulated 
arrays [{t/At)mod M] *A* [k] [{{t-s)/At}mod M] 
in an array for the correlation function (A/j(s)A^) for all 
relevant values of s. At the end of the run the result was 
divided by the number of points accumulated. Multi- 
point and multi-time correlation function are evaluated 
in an analogous fashion. 

Good statistics are difficult to obtain for the higher- 
order correlation functions since the functions are the 
average of a product of multiple factors of the linear den- 
sities For example, the three-point correlation func- 
tions are constructed by averages of quantities which are 
typically on the order of N^, whereas the final average it- 
self is of 0{N). In order to optimize the sampling, many 
relatively short runs of duration R — AM At were per- 
formed and averaged on the fly. The strategy of using 
many short runs seems to be better than the alterna- 
tive of performing a single long run of equal total length, 
perhaps because it reduces the effect of abnormally large 
points which contaminate the signal for a long time. 

Further improvement of the statistics of the calcu- 
lated correlation functions is possible by exploiting the 
isotropy of the system. To simplify the comparison be- 
tween theoretical predictions and the simulation results, 
all wave-vectors were taken to be co-linear along the x 
axis so that k • q = kq, where k — |k| and q = |q|. 



Since the wave-vectors k and q are parallel, the quan- 
tities {Aki:-qi{t)Ags,{t)Al.), {Aky-qy{t)Agy{t)Al.) and 
{Akz-qzit)Aq£{t)Al~) can be computed from the sim- 
ulation in a periodic, cubic simulation box and aver- 
aged to obtain improved statistics. In addition, for 
many of the correlation functions considered here, such as 
{Tk-q{t)Nq{t)T^) , the number of points used to calculate 
the higher-order correlation functions can be effectively 
doubled by averaging over the transverse directions y and 
z. 

The estimation of statistical uncertainty in the sim- 
ulation data is problematic as it involves constructing 
an auto-correlation function for each point measured in 
the time-correlation function p^ . Such a procedure is 
both memory and computationally intensive, and slows 
down the simulation dramatically. In fact, most of the 
computational time of the simulation is spent accumulat- 
ing data and calculating the correlation functions rather 
than performing the molecular dynamics. We therefore 
adopt a simpler approach to estimate the error using the 
symmetry properties of the correlation functions. From 
reflection symmetry, it follows that all correlation func- 
tions are either real or imaginary [ p^ . For a real corre- 
lation function, the imaginary part vanishes and hence 
the imaginary part calculated from the simulation gives 
a rough estimate of the error in the real part, as both 
are calculated from the same configurations and involve 
terms of similar structure. To approximate the statistical 
uncertainty for a real correlation function, a histogram of 
the values of the imaginary part is constructed to deter- 
mine the interval of values containing 96% of the points. 
The size of this interval provides an estimate of the error 
in the real part, taken to be constant for all times of the 
correlation function. Such an approximation seems rea- 
sonable given that the variations in the imaginary part 
in the simulations are observed to be relatively constant 
over the total time interval MAt. For an imaginary cor- 
relation function, the analogous procedure is done using 
the variations in the real part. 

The simulation results for the two-point, two-time cor- 
relation functions were checked against generalized En- 
skog theory results |^l[ . The statistical uncertainty in the 
normalized correlation functions (as obtained by the pro- 
cedure above) are quite small (of the order of 0.001). The 
numerical value for the shear viscosity extracted from the 
exponential decay of the auto-correlation function of the 
transverse velocity Tk was compared to the kinetic theory 
prediction for this quantity pl| ] and excellent agreement 
was observed. 

The time convolution integrals in the mode coupling 
expressions for the higher-order correlation functions 
were evaluated by numerically integrating data for the 
two-point, two-time correlation functions G'^{t) obtained 
directly from the simulation. Since the error bars of the 
G^^lt) are very small, the level of uncertainty in the the- 
oretical prediction for the higher order corelation func- 
tions is negligible in comparison to the uncertainty in 
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the simulation data for the higher-order correlation func- 
tion. Furthermore, no significant differences were noted 
in the convolution integrals calculated using the simula- 
tion data and calculated from high quality functional fits 
of the integrands. In principle, one could also use the 
hydrodynamic forms for all two-point, two-time correla- 
tion functions in combination with an accurate equation 
of state and kinetic theory results for the transport co- 
efficients, but since the simple correlation functions were 
obtained with great accuracy in the simulation, the ac- 
tual data was used. 

As described in the appendix, the dissipative part of 
the vertices for M^^-^, M'^"'-'^ and M'^'^'" have free 
parameters Vh and vth which must be fitted to the 
data if the dissipative contributions are to be included in 
the predictions for both the multi-point and multi-time 
correlation functions. In practice, this is accomplished 
by selecting particular wave-vector magnitudes k and q 
and fitting the parameters according to the simulation 
results. This procedure is illustrated in Fig. 1 for the 
multi-point correlation function 



C 



TL:T 



(t) = {Tk-q{t)Lq{t)T-k)l{N)mkBT 



for wave- vectors k — ko and q = 2fco- Note that al- 
though the dissipative contribution to the overall correla- 
tion function in Eq. ( pO| ) depends on the two parameters 
Vn (N-coupling) and Vh (H-coupling) , these parameters 
can be uniquely determined since the asymptotic time 
behavior is determined entirely by the N-coupling contri- 
bution. Once this parameter is set (here, u„ — —0.18), Vh 
can be determined by fitting the height of the first peak 
(found to be w/j — 0.90). A similar procedure is used for 
the M'^'^'^ vertex in Eq. (17) which is relevant for the 
correlation function 



C 



TT-N 



(t) = {T,_,{t)T,{t)N_,)/S{k), 



and it is found that vth = —0.62. Note that it is, in fact, 
the additional couplings which arise at dissipative order 
that account for the slow decay of the three-point correla- 
tion function in Fig. 1. It is therefore quite apparent that 
ordering of terms using the wave-vector must be done 
carefully for the system under consideration, since con- 
tributions which appear at higher-orders of wave-vector 
can actually dominate lower order terms. 

With the coupling parameters fixed by the fitting pro- 
cedure, one can then compare the simulation results 
with the theoretical predictions for arbitrary wave- vector 
combinations. In Fig. 2, the simulation and theoreti- 
cal predictions for the multi-point correlation functions 
(jTL;T^^-^ and C'^'^'^ [t) are shown as a function of time 
for a number of wave-vector combinations. The remark- 
able agreement between the simulation results and the 
theoretical predictions of both three-point correlation 
functions over all time-regimes and wave-vector combi- 
nations is a clear indication that the formulation of the 
mode-coupling theory is sound. 



It is interesting to see how the theoretical predictions 
of the present formalism compare to those obtained from 
a mode coupling theory in which Gaussian statistical be- 
havior is assumed in the multi-linear basis set. This type 
of assumption roughly corresponds to Kawasaki's origi- 
nal formulation of mode coupling theory [p^ , which is 
based upon a non-linear Langevin equation with Gaus- 
sian noise (or fluctuating forces). Such a Gaussian the- 
ory for the multi-point correlation functions differs from 
the present formulation in two significant ways: First, 
since the subtractions in the multi-linear basis set (see 
Eq. (^) involve static three-point correlation functions 
of the linear densities which vanish under the assump- 
tion of Gaussian statistics, the Gaussian theory neglects 
terms of the form 

(ik(t)i-k) • K^l ■ (i-kik-qiq) 

which appear, for instance, in the mode coupling expres- 
sion for C^^'-^ [t). These terms make an important con- 
tribution to the multi-point correlations functions on all 
time scales, and particularly for short times. Second, 
since the subtraction terms vanish in the Gaussian the- 
ory, the coupling vertices are significantly affected. For 
example, looking at the Eulcr order contributions to the 
vertex M'^'^'^ , the vertex in the Gaussian approximation 
becomes V'^'^''^ , according to 



= ((Tfe:,T,)7V_,) /Sik) EE V^^--^. (26) 

Similar differences between the vertices V in the Gaus- 
sian approximation and the M vertices appear in the 
dissipative parts of the M-^^'-^ and M'^^''^ (see Tables I 
and II). 

In order to assess how each of these differences affects 
three-point correlation functions, we once again consider 
the correlation functions C'^'^'^(t) and C'^^'^{t). The 
first correlation function differs not only in the explicit 
form of the coupling vertices but also in the form of the 
expressions due to the subtraction terms in the basis set. 
The second correlation function, on the other hand, van- 
ishes at i = by symmetry and does not contain con- 
tributions which are directly proportional to two-point, 
two-time correlation functions (see Eq. [T^ ). For this cor- 
relation function, the differences between the Gaussian 
and full mode coupling theory arise solely due to differ- 
ences between the Gaussian vertices V'^^''^ , yTH-.T .^^^^ 
their full counterparts. In Fig. 3, the Gaussian and full 
mode coupling expressions are compared to the simula- 
tion data for the wave- vectors k — ka, q = S/cq- From 
these plots, it is clear that the Gaussian theory poorly 
predicts the time-dependence of C^^''^{t) on all time 
scales, and also gives worse results for the correlation 
function C'^^'^{t). Similar behavior can be seen for other 
wave- vector combinations. 

Turning now to the multi-time correlation functions, 
the simulation results and theoretical predictions for 



7 



the multi-time correlation functions G^^'^ {ti,t2) and 
qTNT fj-^^^^^ for several different wave- vectors are plot- 
ted in Fig. 4 as a function of time t for the time com- 
binations (^1,^2) of {t,t), {t,3t) and {3t,t). The excel- 
lent agreement between the full mode coupling theory 
and simulation results strongly suggests that the assump- 
tions discussed at length in Ref. [|l5[ of what determines 
whether a correlation function decays quickly or not are 
appropriate. These assumptions are necessary to obtain 
mode coupling equations which are local in time. 

Note that the time symmetry properties are evident 
in the two graphs in the first row of Fig. 4, which cor- 
respond to the wave- vectors k — ko, q = 2ko. For these 
wave- vectors, the time symmetries can be obtained by 
noting that 

(Tfe_,(ti -Hi2)A^(ii)r_fe) = (rfe_,(t2)A^r_fe(-ii)) 

since the equilibrium distribution function is stationary. 
Inverting all time arguments and using the properties of 
the densities under time-reversal, one obtains 

(r,_,(i2)A^r-fe(-il)) = 7a(Tfe-,(~t2)A^T_fe(ii)), 

where 7a = 1 for a = iV, _ff and 7a = —1 for T, L. When 
k — q ^ — fc, we find that 

{T^kih + t2)L2k(tl)T^k) = -{T^k{tl+t2)L2k{t2)T^k). 

which implies this correlation function is anti-symmetric 
under interchange of ti with t2 (and therefore vanishes 
when ti = ^2), while {T^kih + t2)L2k{ti)T^k) is sym- 
metric under the exchange of and ^2- It is reassuring, 
though not surprising, that the mode coupling theory re- 
spects these time-reversal properties. 

One may also calculate the multi-time correlation func- 
tions via Eqs. ( p^ ) and (|2^) using the simulation data for 
the multi-point functions G^'^'^ and G^'"'^ . However, 
since the mode coupling results for these functions are 
already in excellent agreement with the simulation data, 
the improvement obtained using the simulation results 
for the G^^ is generally statistically negligible. Further- 
more, the simulations to calculate the multi-point func- 
tions are computationally-intensive compared to calcula- 
tions of the two-point functions. It is therefore far easier 
to generate predictions with small statistical uncertain- 
ties using the mode coupling theory expressions for the 
multi-point functions. 

Since the mode coupling formalism relates the multi- 
time correlation functions to multi-point correlation 
functions, the deficiencies in the Gaussian theory for 
the three-point functions are carried over to the predic- 
tions for three-time functions. This point is confirmed 
by the difference in the behavior of the Gaussian ver- 
sus full mode coupling theory results for the multi-time 
correlation functions shown in Fig. 5. Once again, the 
Gaussian theory predictions for multi-time correlations is 
qualitatively incorrect on all time scales, and particularly 
so for correlation functions which do not vanish when 



ti = t2 = 0. Furthermore, as might be expected from 
the discussion above of the dissipative contribution to the 
three-point correlation functions, the inclusion of the ad- 
ditional couplings arising at dissipative order is essential 
if quantitatively accurate predictions for the multi-time 
correlation functions is desired. 

In principle, in the limit of very small wave-vectors, 
one might expect that the additional couplings in the 
higher-order correlation arising from the dissipative part 
of the vertices become less important and may be ne- 
glected. In fact, this is not always the case since the 
overall order in wave-vector of the various terms in the 
expressions for G^^ is determined by a wave-vector de- 
pendent pre-factor (the vertex) multiplied by the time- 
convolution of two-point, two-time correlation functions. 
The time-convolution of functions such as G'^^(i), which 
vanish as /c — + 0, can give additional factors of the wave- 
vector. Thus, for instance, the contribution from the first 
term in Eq. (|l^) (with a vertex of Euler order) is in fact 
of the same order of magnitude as the contribution from 
the last two terms, which involve vertices of dissipative 
order in the hydrodynamic limit. 

To obtain smaller wave-vectors in a simulation to nu- 
merically check these considerations for dense systems 
in which the mode coupling effects are important, one 
would need to simulate larger systems with more par- 
ticles. There are two difficulties with the simulation 
method applied to larger systems which make it difficult 
to obtain good statistics for the higher-order correlation 
functions. First, since the use of cells is memory intensive 
and the optimal number of cells scales as the cube of the 
length of the system, one must utilize a cell-structure 
for the simulations which is not optimal, leading to a 
reduction in simulation efficiency. Second, the quality 
of the statistics for the higher-order correlation functions 
decreases essentially as the square of the number of parti- 
cles. It is therefore computationally challenging to obtain 
accurate simulation results for the higher-order correla- 
tion functions for larger systems. 



V. SUMMARY AND CONCLUSIONS 

In this paper, the predictions for higher-order corre- 
lation functions based on the mode coupling formalism 
developed in Ref. jl^ were evaluated in the hydrody- 
namic limit for a hard sphere system at moderate densi- 
ties and compared to simulation results. It was demon- 
strated that the mode coupling theory yields excellent 
results for all higher-order correlation functions provided 
that dissipative as well Euler order vertex coupling terms 
are included in the theory. The good agreement between 
the theoretical predictions and the simulation results con- 
firms that the assumptions underlying the mode coupling 
theory of how slow and fast decay of arbitrary densities 
can be separated in a systematic fashion are quite rea- 
sonable. 
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In contrast to some mode coupling theories of simple 
liquids [|l8|,^ 26 2^ , the present mode coupling theory 
includes all multi-linear densities in the set of slow vari- 
ables, does not neglect corrections to the "factorization" 
approximation, and does not assume Gaussian statisti- 
cal properties of the random noise or fluctuating force. 
As the formalism allows exact expressions to be obtained 
for all correlation functions in the thermodynamic limit, 
it provides a systematic way to examine the various as- 
sumptions which must be made in order to predict the 
time-dependence of simple and higher-order correlation 
functions, or to form comparisons with other theories. 
Along these lines, it was demonstrated that the "non- 
Gaussian" behavior of the random noise is important 
for the proper description of the multi-point correlation 
functions on all time scales. In particular, the Gaussian 
theory for three-point functions leads to over-simplified 
coupling vertices which have significant quantitative con- 
sequences and, more importantly, neglects important 
couplings to linear densities. Since the mode coupling 
theory expresses the multi-time correlation functions in 
terms of two-time, higher-order correlation functions, the 
Gaussian theory has similar deficiencies in describing the 
three-time correlation function of linear densities. 

The calculation of higher-order correlation functions 
of extensive linear densities in the hydrodynamic regime 
at low to intermediate densities is computationally in- 
tensive. The poor statistics obtained from the simula- 
tion arises from averaging quantities of order to ob- 
tain a signal of order N. However, since densities of 
tagged particles do not scale with the number of parti- 
cles, higher-order correlation functions of tagged parti- 
cle densities should not suffer from this problem. The 
extension of the mode coupling theory of higher-order 
correlation functions to non-extensive densities of tagged 
particles is straightforward, and will be presented in a 
future publication. 

It is obviously desirable to apply the mode coupling 
formalism to dense and super-cooled liquids where cor- 
relation functions exhibit more complicated time behav- 
ior. In dense systems, there is compelling evidence ||3l| ] 
which suggests that the eigenmode spectrum of the Li- 
ouville operator for simple liquids changes, and a gen- 
eralized "heat" mode becomes long-lived even at fairly 
large wave-vectors. At large wave-vectors, this mode 
roughly corresponds to a self-diffusion mode which 
is slow in dense liquids due to particle caging effects. 
Within the mode coupling formalism, the emergence of 
this short-wavelength collective mode implies that the 
cut-off wave- vector kc for the heat mode becomes on the 
order of inverse molecular length scales. Under these 
circumstances, the mode coupling correction terms to 
the expressions for the higher-order correlation functions 
are not expected to be small and must be considered. 
Appropriately-defined higher-order correlation functions 
may be quite useful in examining the microscopic ori- 
gins of complex relaxation behavior and dynamical het- 
erogeneities. To this end, one may examine the higher- 



order correlation functions at much larger wave-vectors 
using a mode coupling theory in which the modes form- 
ing the basis set for the long-time behavior are associated 
with physical processes on these length scales. In fact, 
the structure of the mode coupling theory suggests that 
measures of dynamical heterogeneity based on multi- 
point correlation functions |lO[ are quite closely related to 
measures based on multi-time correlation functions [ pT| . 
These issues are currently being pursued. 
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APPENDIX A: EVALUATION OF THE 
VERTICES 

In this appendix, all vertices used to formulate numer- 
ical predictions for higher-order correlation functions are 
given for the sake of completeness. To leading order in 
the wave- vectors, all vertices are either of Euler order (or- 
der k), or of dissipative order (order fc^). Since the sec- 
ond term in the expression for the vertices in Eq. (|^) in- 
volves two time derivatives of hydrodynamic densities, it 
is at least of quadratic order in the wave- vectors. There- 
fore the Euler order of any vertex is given by the static 
correlation function (first term) of Eq. (^). This static 
correlation function is imaginary and an odd function of 
wave-vector. The leading order of a vertex of quadratic 
order in wave- vector is therefore given by the second term 
in Eq. (^. For the hard sphere system, all static correla- 
tion functions in the zero wave- vector limit can be evalu- 
ated exactly if the radial distribution function at contact 
g{a) = X is known. The calculation of the vertices at Eu- 
ler order is facilitated by considering the identity, valid 
in the canonical and grand canonical ensembles. 



{Ab)^(3-H{A,B}), 



(Al) 



which links the time derivative to a Poisson bracket of 
the densities. It follows from A = {A , H} and from the 
form of the distribution function: 



dA^dn 

dq dp 



J{A,n}Be-'^'^dr ^ J 

which, by partial integration, yields 
^p-^ j{A,B}e-^'^dT 



dA^dH 
dp dq 



op 



oq \ op 



dV 
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To evaluate the higher-order correlation functions in 
the text, the vertices M^"^'^ and M'^^''^ are needed. The 
latter is the simplest, as Q^^ — Tk-qLq, so 



M, 



TL;T 
k— q,q;k 



-/3-^({Tk-qiq,Tk*})/(Tkrk*) 

- zA;/3-i(rk-qTk*_q)/(rkrk*) = ifcr 



where we used {AB , C} = A{B ,C} + {A, C}B, and the 
facts that {Tq , r*} = and {Lq , T*} = ikT*_^. 

It is straightforward to show that Q^"^ = Tk-q^q — 
|mi?k, so using the above result for {L,T}, we obtain 



M, 



TT-L 
k— q,q;k 



-ik(3- 



{{E^..K}). 



The energy can be split into a kinetic and a potential 
part. The kinetic contribution to M'^'^'^ is easily calcu- 
lated, and M'^'^'^ can be expressed as. 

The second term on the right-hand side of the equation 
above can be evaluated by noting that 

= -i ^ i(k • r„y)x • d,^V{rraj) + 0{k^), 



where r^j — fm ~ , which implies 

{{El°\ LU) = -\zk (N) pjr{x- f)x ■ {d,Vir))g{r) 



dr. 



where g{r) is the radial distribution function. Performing 
the angular integration and writing g{r) = ft,(r)e~'^^^''^ 
so that a partial integration can be performed, leads to 



{{El°\Ll}) 



-ik 



2ttp (N) 
3^3 



dr{r'^h{r))dr = -ikbpNx, 



where b = 27ra'^/3, and x is the radial distribution func- 
tion at contact, x can be estimated using the Carnahan- 
Starling equation of state ||2^ and the expression for the 
pressure p of a hard-sphere system, 



p 



l + rj + rf — rf 



1 



(1 - n? 



where 77 is the packing fraction given by 77 = npa'^ /6. 
Combining all terms, one obtains 

<-'qq;k = *!WP- 

Turning now to the calculation of the dissipative part 
of vertices, their specific wave- vector dependence is deter- 
mined as follows: The derivative of a conserved density 
can be written as 



Ay. — ik - 



where is the current associated with the hydrody- 
namic variable a. What is needed in Eq. (j^) is the dis- 
sipative current = (1 — 'P)J^. Looking first at the 
vertex M'^^''^ , using 



^TH rr. TJ (Tk-q^fqT'-k) rj, rp 

>k_q,q - rk_qi7q - ^^fc^y " T^ 



the vertex can be expressed as 



TH;T 
k-qq 



'kxikx — (Ix)Vh ~ kxQxVfi 

m{N)kBT 



dt. 



(A2) 



where 



Vh 



m{N)kBT 

-(Tk-c^Wj^k) 

~" 



m{N)kBT 



dt 



dt. 



(A3) 



Similar expressions can be obtained for the parameters 
Vn and Vth appearing in the M'^-^''^ and M'^'^'^ vertices 
(see Table II). To obtain the leading behavior for small 
wave- vectors, the wave- vectors in the integrals can be set 
to zero, and the projected dynamics Liouvillian in the 
exponent in Eq. (HQ) can be replaced by the full Liouvil- 
lian. Then, the Green-Kubo expression for the viscosity 
77 can be recognized in the last term of Eq. ([A^), 



V 



{fit)f)dt. 



(A4) 



For the viscosity rj and the heat conduction A (which fig- 
ures in the expression for Al'^'^'^), we take the Enskog 
expressions [ pT| , 

V = Vobp + ^ + 0.7614 6px) , 

A = Xobp ( ^ + + 0.7574 bpx 
V5 bpx 

where the Boltzmann value of the shear viscosity 770 and 
thermal diffusivity Aq are given by 



Ao 



16a2 

75 

64a2 



For the particular parameters of the simulation, it was 
checked by studying the decay of simple correlation func- 
tions that the Enskog expression are accurate. 

In principle, integrals of time-correlations functions of 
products of two currents and a density, as in the expres- 
sion for Vh and v'y^ in Eq. ( [A3| ) can be written in the hy- 
drodynamic limit in terms of transport coefficients and 
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derivatives of transport coefficients with respect to ther- 
modynamic quantities hke the temperature and chemical 
potential. Dissipative contributions such as these have al- 
ready been evaluated by Lim ||3^ in the zero wave-vector 
limit in the context of generalized hydrodynamics. In 
fact, the expression for can be related to the viscosity 
M as 

3 mp' 

whereas Vh can be expressed in terms of the viscosity 
and the derivatives of viscosity with respect to the tem- 
perature and chemical potential. Using the form for u^, 



TH:T 
k— q,q;k 



can be written as 



TH-T 
k— q,q;k 



Vh 



2 7] 

3 mp 



(A5) 



Since it is not known how well the derivatives with re- 
spect to temperature and energy of the approximate ki- 
netic theory expressions for the transport coefficients cor- 
respond to their actual values, Vh, Vn and Vth are taken 
as free parameters which will be fitted from simulation 
data. 

The expressions for the vertices that are needed in the 
text are listed in Tables I and II. Also tabulated are 
the vertices one would obtain from a Gaussian theory 
in which static three point correlation functions are set 
to zero. These Gaussian vertices are denoted by V21. 
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Figure Captions 



Figure 1. The fitting procedure for the three-point correlation function C"^^'^(t) for the wave- vectors k = kg, q = 2ko, 
where koa = 0.398867. The un-connected circles are the simulation results, the solid line is the full mode coupling 
results, the dotted line is the mode coupling results with Euler vertices, and the long-dashed and dot-dashed lines 
represent the contributions from the iV-dissipative and if-dissipative vertices, respectively. For clarity, the statistical 

uncertainties in all quantities have been omitted. 

Figure 2. The correlation functions C"^"^'^ {t) (left panels) and C^^'''^{t) (right panels) as a function of reduced time 
at various wave-vectors. In the top row, the wave-vector arguments are k = ko, q = 2fco (open un-connected circles: 
simulation results, solid line: MCT prediction), and k = 2kQ, q = k^ (open im-connccted squares: simulation results, 
dotted line: MCT prediction). In the middle row, the wave-vector arguments are k = ko, q = S/cq (open un-connected 
circles: simulation results, solid line: MCT prediction), and k = 3fco, q = ko (open un-connected squares: simulation 
results, dotted line: MCT prediction). In the bottom row, the wave- vector arguments arc k = 2fco, q = 3fco (open 
un-connected circles: simulation results, solid line: MCT prediction), and k = 3ko, q = 2ko (open un-connected 
squares: simulation results, dotted line: MCT prediction). For clarity, the statistical uncertainties in all quantities 
have been omitted. 

Figure 3. The full mode coupling theory (MCT), Gaussian MCT, Euler order MCT predictions and simulation data 
for the correlation functions C'^'^'^ [t) (left panel) and C'^^'^{t) (right panel) at wave-vectors k = ko and q = 3ko- 
In both panels, the un-connected circles are the simulation data, the solid, dotted and dashed lines represent the full 
MCT theory, the Euler order MCT theory, and the Gaussian MCT theory results, respectively. The error estimates 
represent 96% confidence intervals. Note that the Gaussian MCT theory is qualitatively incorrect on all time scales 
for (7^^;^(t). 

Figure 4. The multi-time correlation fimctions G'^^'^{ti, t2) and G"^^"^(ti, t2) as a fmiction of reduced time for various 
wave- vector combinations. In all panels, the un-connected dots, x 's and triangles correspond to the simulation results 
for the time arguments ti = t, t2 = t, ti = 3t, t2 = t and ti = t, t2 = 3t, respectively. The solid, dashed, and dotted 
lines correspond to the respective mode coupling predictions. The results in the top, middle and bottom rows arc 
for the wave- vector arguments k = ko, q = 2ko, k = ko, q = 3fco, and k = 2ko, q = ko- For clarity, the statistical 
uncertainties in all quantities have been omitted. 

Figure 5. The full mode coupling theory (MCT), Euler order MCT, and Gaussian MCT predictions and simulation 

data for the multi-time correlation fimctions G'^^'^ (ti.t2) (left panel) and G'^ ^'^ (t i . t2) (right panel) at wave- vectors 
k ^ ko and q — 3fco and time arguments ti — 3t, t2 = t. In both panels, the un-connected circles are the simulation 
data, the solid, dotted and dashed lines represent the full MCT theory, the Euler MCT theory, and the Gaussian 
MCT theory results, respectively. The error estimates represent 96% confidence intervals. 
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TABLE I. Expressions for the leading behavior of the Euler vertices M and their Gaussian counterparts V' 
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TABLE II. Expressions for leading behavior of the dissipative vertices M^^ and their Gaussian counterparts V'^^ . Note that 
in the table k and q stand for the x component of k and q, respectively. 
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